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An  approximate  solution  is  presented  for  the  spherical  diffusion  equation  for  the  spherical  particles  in  a 
physics-based  lithium-ion  battery  model.  This  approximate  solution  is  compared  to  different  numerical 
and  analytic  solutions  for  various  boundary  conditions.  These  comparisons  reveal  that  our  approximate 
solution  can  provide  accuracy  and  time-efficiency  in  simulation.  This  approximate  solution  is  much  faster 
than  the  numerical  and  truncated  analytical  solutions  at  high  current  rate,  and  shows  better  long-time 
accuracy  than  the  short-time  analytical  solution.  This  approximate  solution  can  also  be  used  as  the  porous 
electrode  model. 
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1.  Introduction 

Physics-based  models  are  widely  used  to  predict  the  behav¬ 
ior  of  lithium  ion  cells  and  batteries  [1].  In  a  typical  cell  model, 
each  porous  electrode  is  described  by  a  pseudo  two-dimensional 
domain  with  solid-phase  diffusion  occurring  in  spherical  particles, 
for  example,  at  each  local  position  in  a  porous  electrode  (see  Fig.  1  of 
Ref.  [1]).  This  solid  state  diffusion  occurs  because  lithium  ions  are 
intercalated  or  deintercalated  through  electrochemical  reactions 
at  the  surface  of  the  solid  phase  particles  [2,3].  For  this  reason,  the 
concentration  of  lithium  ions  at  the  particle  surface  is  an  impor¬ 
tant  variable  which  is  needed  in  the  simulation  to  determine  the 
reaction  rate,  for  example. 

Several  numerical  methods  have  been  used  to  solve  the  solid- 
phase  diffusion  equation  in  the  r  dimension  of  spherical  particles 
(the  pseudo  second  dimension)  and  obtain  the  surface  concen¬ 
tration  of  lithium  ions.  For  example,  the  finite  difference  (FD)  [4] 
and  orthogonal  collocation  on  finite  elements  (OCFE)  [5,6]  meth¬ 
ods  have  been  used.  These  methods  require  discretization  of  the 
particle  radius  into  node  points  or  elements  with  interior  points; 
and,  consequently,  one  or  more  ordinary  differential  equations 
(ODEs)  are  obtained  at  each  point  or  element  (this  method  is  often 
called  the  method  of  lines).  Flowever,  the  computational  error 
caused  by  the  discretization  of  a  continuous  domain  increases  dra¬ 
matically  if  there  are  large  concentration  gradients  at  the  surface 
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of  the  particles.  If  too  many  node  points  or  elements  are  applied  to 
reduce  the  error  in  discretization,  the  simulation  is  severely  slowed 
due  to  the  large  number  of  equations. 

The  analytic  solution  to  the  spherical  diffusion  equation  includes 
an  infinite  series.  To  evaluate  this  infinite  series,  the  series  is 
approximated  by  truncating.  The  truncated  series  includes  a  finite 
number  of  terms.  Evaluation  of  these  terms  requires  significant 
computation  time  for  cases  with  large  concentration  gradients  at 
the  surfaces  of  the  particles,  which  can  occur  during  pulsing  for 
example. 

An  approximate  analytic  solution  is  presented  in  this  paper 
to  minimize  the  truncation  error  associated  with  implementation 
of  the  series  with  a  finite  number  of  terms  included  in  the  ana¬ 
lytic  solution.  It  is  shown  that  only  a  few  terms  are  needed  in  the 
approximate  solution  together  with  an  additional  term  compen¬ 
sating  for  the  truncation  error.  The  resulting  approximate  solution 
yields  accurate  results  with  less  computation  time  than  required 
for  numerical  methods  or  a  large  number  of  terms  in  the  analytic 
solution. 


2.  Model  development 


The  diffusion  of  lithium  ions  in  a  spherical  solid  phase  particle 
follows  Fields  law  and  is  described  by  the  following  partial  differ¬ 
ential  equation  (PDE)  in  the  spherical  coordinates: 


dc 

dt 


D 12. 
r2  or 


(1) 


at  t  =  0  for  0  <  r  <  R  c(t  =  0)  =  Co 


(2) 
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Nomenclature 

List  of  symbols 

c  concentration  of  lithium  ions  in  the  solid  particles 

(molm-3) 

Cmax  maximum  concentration  in  particle  (mol  m-3 ) 

C  dimensionless  concentration  of  lithium  ions 

Cs  dimensionless  concentration  of  lithium  ions  at  par¬ 

ticle  surface 

C  dimensionless  average  concentration  of  lithium 

ions 

D  solid  phase  diffusion  coefficient  for  lithium  ions  in 

the  particles  (m2  s-1) 
ejv  truncation  error 

j(t)  molar  flux  of  lithium  at  the  particle  surface 

(molm-2  s-1) 

n  index  of  eigenvalues  and  eigenfunctions 

N  number  of  eigenfunctions  kept  in  the  series 

Qn  eigenfunction 

Q,?px  approximated  eigenfunction 

r  radial  coordinate  (m) 

R  particle  radius  (m) 

r  dimensionless  radial  coordinate 

t  time  (s) 

8{  r)  dimensionless  molar  flux  of  lithium  ions  at  the  par¬ 
ticle  surface  in  the  f  direction 
Xn  eigenvalue 

r  dimensionless  time 


dc 

dr 


=  0 

r= 0 


=Kt) 

r=R 


(3) 

(4) 


where  c  is  concentration  of  lithium  in  the  solid  particles,  D  is  the 
solid  phase  diffusion  coefficient  for  lithium  in  the  particles,  j(t)  is 
the  time  dependent  boundary  flux,  and  R  is  the  radius  of  particle. 
Assuming  that  the  diffusivity  D  and  the  particle  radius  R  are  con¬ 
stants,  these  model  equations  can  be  rewritten  in  the  dimensionless 
form  by  choosing  the  following  dimensionless  variables: 


(5) 

Cmax 

Dt 

(6) 

X=R2 

_  j(t)R 

D(t)cmax 

(7) 

r 

r  ~  R 

(8) 

where  cmax  is  the  maximum  lithium  concentration  in  the  particle. 
The  dimensional  equations  (1)  through  (4)  can  be  converted  into 
the  following  dimensionless  equations  by  using  Eqs.  (5)-(8) 


dC  _  J_d_  (-2dC\ 
9r  “  Pd?  \  df) 


(9) 


C(r  =  0)  =  C0 


9C 

dr 


=  0 

r= 0 


dc 

d? 


r=  1 


-8(r) 


(10) 

(11) 

(12) 


The  analytic  solution  to  these  equations  is  [7] 
C(f,  r)  .  Co  +  ^  J(r)  -  ^P2 


-3  + 

Jo  Sm(An) 

_2ysmO;nf)  -k2T  fTS{r,)exir'dr,  (13) 

sin(X„)  J0 


The  last  two  terms  in  Eq.  (13)  can  be  combined  by  letting 
Qn(r)  =  -2e~x"r 
Eq.  (13)  then  becomes 


A"r  [  <5(r')e*"T'dr' 

Jo 


C(r,  r)  =  C0  +  T5(T)-^r2-3 


[  S(r')dr' 

Jo 


(14) 


Esin 

r  si 

n= 1 


sin(A.nr) 
r  sin(An) 


Qn('t)  + 


28(r) 

A2 


(15) 


Eq.  (14)  can  be  written  as  a  differential  equation  by  taking  the 
derivative  of  each  side  of  Eq.  (14)  with  respect  to  t: 

=  -^nQn(t)  -  28(t)  for  n  =  1  to  oo  (16) 

with 


Qn(r  =  0)  =  0  for n  =  1  to  oc  (17) 

where  An  is  the  nth  eigenvalue  calculated  from  the  following  equa¬ 
tion 


Xn  —  tan Xn  =  0  n  =  l,2,  ••• 

Note  that  when  <$(r)  is  a  constant,  Eq.  (16)  yields 

Qn(r)  =  -^[1  -  exp(-X2r)] 

An 

at  the  surface  of  the  particle 
Cs(r)  =  C(f=l,  r) 


(18) 


(19) 


(20) 


Eq.  (15)  can  be  used  to  write  the  dimensionless  surface  con¬ 
centration,  Cs(t),  in  terms  of  the  average  concentration  and  the 
eigenfunction,  Q„: 


Cs(r)  =  C(r)  +  ^Q„(r)  (21) 

n= 1 

where  C(r)  =  3  J1  r2C(r,  r)dr  is  the  average  concentration  and  is 
determined  by 

^  =  -3S(r)  C(r  =  0)  =  C0  (22) 

Eq.  (21)  presents  the  theoretical  form  for  the  complete  solution 
for  Cs(r)  if  an  inhnite  number  of  terms  are  included  in  the  series. 
However,  to  calculate  a  value  for  the  surface  concentration,  one 
normally  retains  N  terms  in  the  series  by  using  a  truncated  solution: 

N 

Cstrun(r)  =  C(r)  +  ^Qn(r)  +  eN(r)  (23) 

n= 1 
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where  the  truncation  error,  eN(r),  between  the  truncated  solution 
given  in  Eq.  (23)  and  the  complete  solution  given  in  Eq.  (21 )  is  given 
by 


ew(t)  =  E  QnM 


(24) 


n=N+ 1 


Now,  let’s  develop  an  approximation  for  this  truncation  error.  To 
do  so,  let’s  find  an  approximate  eigenfunction  Q^px.  That  is,  let 

-vapx  _  28(t)  r 


QnP  =  -- 


xl 


-[1  -  exp(-A„r)] 


(25) 


(note  that  this  is  the  eigenfunction  for  a  constant  8  given  in  Eq.  (1 9)). 
Next,  take  the  derivative  of  both  sides  of  Eq.  (25)  with  respect  to  r 
to  obtain  the  following  equation 

d(2naPX  =  25(f)- 


dr 


X2n  dr 
-  exp(-A^r)]  with  Q^px 


T=0 


=  0 


(26) 


For  large  values  of  n,  the  nth  eigenvalue  Xn  will  be  large  enough  to 
satisfy  the  relationship 


25(f) »  A^)[i  _  exp(-A.^r)] 

and  Eq.  (26)  can  be  simplified  to 

dQ^px 


dr 


-XtQZpx-2S(r) 


Next,  let  the  approximate  truncation  error  be 

oo 

4pxM=  E  Q"apx(r) 


(27) 


(28) 


(29) 


n=N+ 1 


where  N  is  large  enough  to  insure  that  Eq.  (27)  is  valid.  Since  Q^px 
is  a  function  of  n,  the  summation  5^^N+1Q^px  can  converted 
approximately  to  the  integration  of  Q^px  with  n: 


x  pOO 

V  QnaPX  =  /  < 

-N+ 1  Jn+1 


Q„  dn 


n=N+ 1 


For  large  n,  it  is  known  that 
X„  **  (j  +n)  71 

Consequently, 
dXn  =  7i  dn 

Substitution  of  Eqs.  (25)  and  (32)  into  Eq.  (30)  yields 
E  Qnapx  =  -|  r  -  exp(-A2r)]dXn 

n=JV+l  •'AIV+1  n 


(30) 

(31) 

(32) 

(33) 


The  right  hand  side  of  Eq.  (33)  can  be  integrated  analytically  to 
obtain: 


.1  f 

71  h 

J  Am 


28{r) 


[1  -  exp{-X2t)]dXn 


=  -2S(t) 
Since 


l-exp(— X2  jT)  fx  ... 

- 3T - ± —  +  \  -  erfc(Ajv+i  Vx ) 

7lAjy+^  V  71 


f*OQ  OO 

1  _  i  1  ..  ^  1 

^N+1  njx 2  n~  2^  x2 

"'Xn+'  "  n=N+l  " 


(34) 


(35) 


substitute  Eq.  (35)  into  (34)  and  to  obtain 

2S(r)-[l  -  exp(-X2nr)]dXn  =  -25(f)  J  (  E  i  )  [1 


1  f 

71  Aw 


71  /l  A.2 

7  AN+1  n 


\n=N+ 1 


-exp(-X^+1r)]  +  \/  - 


(36) 


The  derivation  from  Eqs.  (34)  to  (36)  effectively  removes  the 
error  generated  from  converting  addition  to  integration  (Eq.  (30)) 
because  Eq.  (36)  satisfies  the  long-time  accuracy.  That  is,  as  oo, 
exp(-A„r)  ->  0,  and  according  to  Eq.  (25), 


napx  _  2 8(r) 

Qn  —  z~2  O.S  r  —x  OQ 

An 

therefore 


oo  oo 

Eo,?px=-2S(t)E  A 

An 


n=N+ 1 


n=N+ 1 


(37) 


(38) 


Also  as  f->-  oo,  y^r/tr)  erfc(Aw+i  Vf )  ->•  0,  and  the  right-hand-side 
of  Eq.  (36)  becomes 

-25(f)  |  (  E  72  J  I1”  exP(-4+1T)]  +  \/A 


Vn=N+l 


oo 

=-25(t)E  i 


n=N+ 1 


(39) 


Comparing  Eq.  (38)  and  (39),  we  obtain  that  the  following  approx¬ 
imation 

oo  (  /  oo  \ 

E  Qi?px  ^  -25(r)  <  (  e  ^  n1  ■  exp(-AN+it)i 


n=N+l 


Vn=N+l 


which  has  good  accuracy  for  oo. 
Since 


as  r  ->  oo 


V-  =  — 

tr*  10 

we  can  write 


oo  oo  N  N 

El  1  1  1  1 

2^x2  _  10  X^X2 

n=JV+l  "  n=l  n  n=l  "  n=l  " 

Substitution  of  Eq.  (42)  into  (40)  yields 


(40) 


(41) 


(42) 


n=N+l 


-exp(-A2  ,t)]  +  \  J- 


(43) 


Thus,  the  approximate  truncation  error  becomes 

4PX(T)  =  ~2S(T)  j  (  ^  “  E^2  j  [1  “  eXP(-XN+lr)l 


+\  -  erfc(AN+1  Vr)  > 
V  JT  J 


(44) 
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Using  Eq.  (44)  to  approximate  the  truncation  error  in  Eq.  (23),  we 
obtain  our  approximate  solution  for  the  surface  concentration: 

N  (  /  N  \ 

Csapx(r)  =  C(r)  +  ^Qn(r)-25(r)<|  (  Jq~Yj2  )  t1 


n= 1 


-exp(-A^  r)]  +  J  -  erfc(AN+ix/r) 


(45) 


when  <$(r)  is  a  constant,  Eq.  (45)  becomes 
Csapx(r)  =  C0-3«r-2«|(^-E^J  [I"  exP(-AN+i^)l 


+a/J 


-  -  exp(— A^r)]  (46) 


when  5(r)  =  -  sin(r)  Eq.  (45)  becomes 

Capx(r)  =  C0-3[cos(r)-l]  +  2sin(r)|  (jq~Yt2  )  I1 


n= 1 


-  exp(-^  T)]  +  J-  erfc(XN+1  Vr ) 


N  r 

+E 

n= 1  L 


2e_;inT  -  2  cos(r)  +  sin(r) 


+ 1 


(47) 


Eq.  (45)  is  our  approximate  solution  for  the  dimensionless  surface 
concentration  as  a  function  of  time.  It  will  be  shown  below  that  our 
approximate  solution  provides  an  accurate  value  for  the  surface 
concentration  with  only  a  few  terms  (N  =  5,  e.g.)  in  the  series  in  Eq. 
(45)  in  comparison  to  the  truncated  solution  (Eq.  (23))  with  a  large 
number  of  terms.  It  is  also  shown  that  our  approximate  solution 
with  only  a  few  terms  agrees  well  with  numerical  solutions  for 
both  constant  boundary  flux  and  time-dependent  boundary  flux 
expressions. 


3.  Results  and  discussion 

Our  approximate  solution  for  the  surface  concentration,  Cfpx(t), 
can  be  obtained  in  analytic  form  when  8  is  a  constant  or  a  simple 
function  of  dimensionless  time  or  by  numerically  solving  N  + 1  ODEs 
(one  for  C  and  N  for  Ch  through  Qn).  The  simulations  presented 
below  were  based  on  the  discharge  of  the  positive  electrode  in  a 
lithium  ion  pouch  cell  (the  design  parameters  for  this  type  of  cell  are 
available  in  Ref.  [8]  and  listed  in  Table  1 ).  Where  the  dimensionless 
flux  equivalent  to  1  C  rate  is  calculates  as 

s'c-bb£;-0M0452  1481 

The  initial  condition  is  C0  =  0.5  and  the  simulation  stops  when  Cs 
reaches  0.95. 


3.1.  Comparisons  with  truncated  solutions 

Comparisons  between  our  approximate  solution  and  the  trun¬ 
cated  solutions  are  presented  in  Fig.  1.  In  this  case,  the  boundary 
flux  8  =  -  20  is  equivalent  to  1 00  C  rate.  As  shown  in  Fig.  1 ,  with  only 
5  terms,  our  approximate  solution  provides  the  same  accuracy  as 
the  truncated  solution  with  5000  terms.  Note  that  the  300  term 
truncated  solution  is  not  accurate. 


Fig.  1.  The  comparison  between  the  reformulated  eigenfunction  solution  and  trun¬ 
cated  solutions. 


3.2.  Comparisons  with  numerical  solutions 

To  compare  our  approximate  solution  to  numerical  solutions, 
two  commonly  used  numerical  approaches  were  used:  the  finite 
difference  (FD)  method  and  the  orthogonal  collocation  on  finite 
element  (OCFE)  method.  As  shown  in  Fig.  2(a),  the  five-term 
approximate  solution  agrees  well  with  FD  solution  1  with  2000 
node  points.  If  the  number  of  node  points  is  reduced  to  100,  the  FD 
solution  2  shows  poor  accuracy.  The  FD  solution  1  deviates  slightly 
from  the  reformulated  eigenfunction  solution  around  the  initial 


b  0.95 
0.9 
0.85 
0.8 
0.75 
0.7 
0.65 
0.6 
0.55 
0.5 

0  0.0001  0.0002  0.0003  0.0004 

7 


—  Approximate  Solution  (N=5) 
o  OCFE  Solution  1  (50  Elements,  2  Interior  Points) 
x  OCFE  Solution  2  (10  Elements,  2  Interior  Points) 
a  OCFE  Solution  3  (100  Elements,  1  Interior  Point) 
x  OCFE  Solution  4  (20  Elements,  1  Interior  Point) 


Fig.  2.  The  comparison  between  the  reformulated  eigenfunction  solution  and  (a)  FD 
solutions,  (b)  OCFE  solutions. 
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Table  1 

The  design  parameters  for  positive  electrode  of  MSA  pouch  cell. 


Symbol 

Description 

Value 

Unit 

he 

1  C  rate  current 

1.656 

A 

Da 

Diffusion  coefficient  in  the  particle 

1.213  x  10“14 

m2  s_1 

Sp 

Total  electroactive  surface  area  of  positive  electrode 

1.167 

m2 

Cmax 

Maximum  concentration  in  particle 

51,410 

mol  m-3 

F 

Faraday’s  constant 

96,487 

Cmol-1 

Rp 

radius  of  particle 

8.5  x  10“6 

m 

a  Evaluated  at  30  °C  using  Arrhenius’s  correlation  with  activation  energy  of  29  kj  mol-1 . 


Fig.  3.  The  surface  concentration  with  periodic  boundary  flux  (a)  low  frequency,  (b) 
high  frequency. 

time  r  =  0,  even  though  these  two  solutions  agree  perfectly  with 
each  other  over  the  most  of  the  time  range.  In  FD  solution  1,  the 
calculated  initial  value  is  Cs  =  0.5067  at  t  =  0,  which  is  not  consis¬ 
tent  with  the  set  value  0.5.  The  reason  is  that,  in  the  FD  method,  the 
dependent  variables  at  the  two  boundary  points  are  determined 
by  algebraic  equations,  and  the  error  is  caused  in  the  initialization 
of  DAE  system.  For  our  approximate  solution,  all  dependent  vari¬ 
ables  are  determined  by  ODEs  and  no  initialization  is  needed,  so 
the  solution  can  start  exactly  at  the  set  value. 

Fig.  2(b)  presents  a  comparison  between  our  approximate 
solution  and  the  OCFE  solutions  with  different  settings.  As 

Table  2 

The  simulation  time  for  different  approaches. 


Approach  Solution  time  (s) 


Approximate  solution  (N  =  5)  (Eq.  (45))  0.015 

FD  solution  1  (2000  points)  4.98 

OCFE  solution  1  (50  elements  and  interior  points)  0.391 

OCFE  solution  3  (100  elements  and  interior  point)  0.406 


shown  in  the  figure,  the  five-term  approximate  solution  agrees 
with  the  OCFE  solutions  1  and  3  in  which  finer  element 
mesh  size  is  applied  and  more  equations  are  included  to 
ensure  sufficient  accuracy.  The  OCFE  solutions  2  and  4  with 
larger  element  mesh  size  and  less  equations  have  significant 
error.  In  OCFE  method,  all  equations  are  ODEs  and  no  ini¬ 
tialization  is  involved,  so  consistent  initial  conditions  can  be 
obtained. 

The  solution  times  for  different  approaches  are  listed  in  Table  2. 
The  reformulated  eigenfunction  solution  is  much  faster  than  the 
other  numerical  solutions  with  the  same  level  of  accuracy.  Being 
time-efficient  is  a  great  advantage  for  the  reformulated  solution 
over  the  numerical  solutions.  It’s  also  found  in  Table  2,  the  OCFE 
method  is  faster  than  the  FD  method,  and  furthermore,  there 
is  no  initial  error  for  OCFE  method;  therefore  in  the  following 
case  studies,  the  OCFE  solution  with  50  elements  and  2  inte¬ 
rior  points  is  used  to  validate  the  accuracy  of  our  approximate 
solution. 


Fig.  4.  The  comparison  between  the  reformulated  eigenfunction  solution  and  short- 
time  analytic  solution  (a)  large  boundary  flux,  (b)  small  boundary  flux. 
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3.3.  Comparisons  with  the  short-time  analytic  solutions 

Using  the  method  described  by  Atlung,  et  al.  [9]  if  the  bound¬ 
ary  flux  <5  is  constant,  a  short-time  analytic  solution  for  Cs  can  be 
obtained  (see  Appendix  A)  as  follows: 

Cs(t)  =  C0-S[  1  -  er  erfc(-Vt)]  (49) 

Fig.  3(a)  shows  a  comparison  between  our  approximate  solu¬ 
tion  and  the  short-time  solution  with  a  large  boundary  flux  S  =  -  20. 
According  to  this  figure,  the  predictions  from  our  approximate  solu¬ 
tion,  the  short-time  analytic  solution,  and  the  OCFE  solution  are 
consistent.  Fig.  3(b)  presents  the  comparison  with  small  bound¬ 
ary  flux  S  =  -  0.2,  which  corresponds  to  a  1  C  discharge  rate.  In  this 
case,  since  it  takes  longer  time  for  Cs  to  charge  from  0.5  to  0.95,  the 
short-time  analytic  solution  loses  its  accuracy  and  deviates  from 
the  other  two  curves  at  the  final  stage,  but  our  approximate  solu¬ 
tion  still  agrees  well  with  the  OCFE  solution  (which  can  be  safely 
considered  as  accurate)  over  the  entire  time  range. 

The  above  comparisons  show  that  the  use  short-time  analytic 
solution  is  only  valid  over  a  small  time  range.  Another  limitation  for 
the  short-time  analytic  solution  is  that  it  only  works  for  a  constant 
boundary  flux.  Therefore,  our  approximate  solution  is  more  useful. 


3.4.  With  time-dependent  boundary  flux 

In  this  paper,  two  time-dependent  boundary  flux  con¬ 
ditions  are  considered:  <$(r)  =  -  [1  +  sin(100r)]  and  <$(r)  = 
-  [1  +  sin(lOOOr)].  The  results  are  presented  in  Fig.  4(a)  and  (b). 
In  Fig.  4(a),  at  a  lower  frequency,  the  five-term  eigenfunction 
solution  shows  good  accuracy  as  validated  by  the  OCFE  solution.  In 
Fig.  3(b),  however,  the  five-term  eigenfunction  solution  becomes 
less  accurate  when  the  frequency  is  increased  to  1000,  and  nine 
eigenfunction  terms  are  needed  for  this  case  (N= 9).  The  reason 
for  the  additional  terms  can  be  explained  by  considering  Eq.  (27). 
Since  the  time-derivative  of  the  boundary  flux  (dS/dr)  increases 
with  frequency,  a  larger  Xn  is  required  to  keep  the  right  hand  side 
of  Eq.  (27)  small  enough;  and,  therefore,  more  eigenfunctions  are 
needed  in  the  series. 

3.5.  Application  in  porous  electrode  model 

Our  approximate  solution  can  also  be  used  to  simply  the  porous 
electrode  model  (the  Pseudo-2D  model)  [1,2].  The  2C  and  20  C 
discharge  profiles  simulated  by  pseudo-2D  model  with  different 
solution  approaches  for  the  solid-phase  diffusion  equations  are  pre¬ 
sented  in  Fig.  5(a)  and  (b),  and  in  each  plot,  the  OCFE  solution  with 
50  elements  and  1  interior  point  is  used  as  the  complete  solution. 
As  shown  in  Fig.  5(a),  at  the  2  C  rate,  both  the  1-term  approximate 
solution  and  the  5-element  OCFE  solution  agree  with  the  complete 
solution  but  the  50-term  truncate  solution  loses  accuracy  at  the 
end  of  discharge.  As  shown  in  Fig.  5(b),  when  the  discharge  rate 
increases  to  20  C,  the  1-term  approximate  solution  still  shows  good 
accuracy  but  the  5-element  OCFE  solution  and  the  50-term  truncate 
solution  deviate  significantly  from  the  complete  solution. 

4.  Conclusion 

An  approximate  solution  is  presented  for  the  surface  concen¬ 
tration  in  a  spherical  particle  with  a  constant  or  a  time-dependent 
flux  boundary  condition  at  the  surface.  Our  approximate  solution 
for  the  surface  concentration  (Eq.  (45))  applies  for  various  bound¬ 
ary  conditions:  large  constant  flux  for  short  time,  small  constant 
flux  for  long  time,  low  frequency  periodic  boundary  flux,  and  high 
frequency  periodic  boundary  flux.  This  approximate  solution  has 
proven  to  be  advantageous  over  the  commonly  used  solutions: 

1 )  It  contains  much  fewer  terms  than  the  truncate  solution  as  the 
truncation  error  is  effectively  approximated. 

2)  It  shows  the  same  level  of  accuracy  with  numerical  solutions 
including  large  numbers  of  node  points  or  mesh  elements,  but 
is  much  faster. 

3)  It  shows  much  better  long-time  accuracy  than  the  short-time 
analytical  solution. 

4)  This  approximate  solution  also  works  for  the  porous-electrode 
model  in  which  the  particle  surface  flux  is  dynamic. 

Therefore,  this  approximate  solution  method  greatly  improves 
the  simulation  efficiency  and  accuracy  of  physics-based  Li-ion  cell 
models. 

Appendix  A.  The  short-time  solution  for  particle  diffusion 

Dimensionless  governing  equation  for  solid-phase  diffusion: 

d£  =  11  (-2  d£\ 

dr  f2  dr  y  dr ) 

Initial  condition: 

C(r  =  0)  =  C0 


(A-l) 

(A-2) 
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Boundary  conditions: 


dc 

w 


=  0 

r= 0 


dc 

W 


=  8 


r=l 


Taking  Laplace-transformation  C(r,  r) 
equation 


(A-3) 

(A-4) 

u(s,  r)  for  governing 


Solving  Eq.  (A-5)  to  obtain 
„  ^  sinh(Vsr)  cosh(Vsr)  ,  C0 

u  =  L2 - _ - b  Li - - - I - 


(A-5) 


(A-6) 


where  C|  and  C2  are  integral  constants  to  be  determined  by  bound¬ 
ary  conditions. 

As  u  cannot  be  infinitely  large  at  r  =  0,  according  to  Eq.  (A-6) 


C,  =0 


(A-7) 


and 


sinh(Vsr)  C0 

z  L2 - - - 1 - 

r  s 

Applying  boundary  condition  at  x  =  1 
=  -C2  sinhtVs)  +  C2Vscosh(Vs) 

r=l 

=  = _ 5 _ 

s  s  [sinh(Vs)  -  Vscosh(Vs)] 


dc 

W 


(A-8) 


(A-9) 


sinh  (Vsr)5  +  C0 

rs  [sinh(Vs)  -  Vscosh(Vs)]  s 


(A-10) 


According  to  Eq.  (A-10),  the  surface  concentration  at  r  =  1  can 
be  rewritten  as 
5 


Us  = _ +  Co 

s  sVscoth(Vs)  -  s  s 

As  s ->  00,  coth(Vs)  ->  1,  and 
us  = 


S  |  Cp 
sVs  -  s  s 


(A-ll) 


(A-12) 


Take  inverse  Laplace  transformation  us(s )  ->  Cs(r)  for  Eq.  (A-12) 
and  obtain 


Cs  =  C0  -  5[1  -  exerfc{-Vr)] 


(A-13) 
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